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SUMMARY 


This report describes several improvements in "state-of-the-art" software 
which were made during the development of the RAINPAL* concept and in subse- 
quent extensions of the concept. One of the improvements is to formulate the navi- 
gation equations in Earth-fixed coordinates, with the coordinate center and axes 
chosen so that . ■ .. /■' . : 

1. minimal auxiliary calculations are required for pilot displays and 

2. high numerical accuracy is retained with single-precision, fixed- 
point arithmetic. 

Use has been made of area navigation (RNAV) waypoints, in addition to runway 
reference points, as coordinate centers during the various phases of flight. 

Another improvement is in the formulation of the variational equations which 
are, used in the navaid data processing algorithm (a modified Kalman filter with a 
square-root implementation). Simplifications in these equations result with the use 
of a new concept which mathematically ties the variational . equati on reference to the 
platform. An added advantage of this new approach is that accelerometer measure- 
ments do not occur in the Jacobian matrix of the' acceleration with respect to the 
state. , . \ ! ' 

Detailed discussion of the specifics of the RAINPAL software are presented 
as well as the extensions of the RAINRAL concept to other phases of flight. 


* RAINPAL is an acronym for "recursive-aided-inertial-navigation-for-precision- 
approach-and-landing." The RAINPAL system concept has been validated through 
flight tests in a recent NASA Ames inhouse and contractual effort. 
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I. INTRODUCTION 

An attractive concept for an all-weather aircraft navigation system 

which can potentially satisfy the requirements for all phases of flight is 

to use measurements from external navigation aids in conjunction with an 

inertial navigation system (INS). The. excellent high-accuracy, short-term 

characteristics of the INS provide a reference for filtering the noise in the 

navigation aid measurements. The position information from appropriate 

navigation aids provides a good reference for removing the drift characteristics 

of the INS. Current "state of the art" hardware for INS, navigation aids 

and onboard computers is potentially applicable for such an aided inertial 

navigation system. The feasibility of this concept for approach and landing 

operations has been demonstrated in' recent NASA/ Ames flight tests of an 

* • 

experimental aided inertial navigation system called RAINPAL. This system 

was developed in a joint NASA/Ames and contractual effort aimed at the " 

definition, formulation and validation of new or advanced navigation techniques 

which are compatible with the computational capabilities of "off the shelf" 

onboard computers. 

Kalman filter theory provided the basis for the onboard computer 

algorithms used in combining the INS and navigation aid data. The square-root 

1,2 

formulation of this filter which incorporates random forcing functions in 
the INS error model forms part of the RAINPAL software. Other advanced 
techniques tested in the RAINPAL software included the solution of the 
navigation equations in a runway- referenced (Earth-fixed) coordinate frame 
and a new formulation of the error state vector for inertial navigation systems. 

* RAINPAL is an acronym for "recursive-aided-inertial-navigation-for- 
precision-approach-and-landing." 
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3 4 5 6 

Previous reports ' ’ ’ have described both simulation results and flight 
test results associated with the overall effort. The mathematical formulations and 
software details of the experimental system have not been presented in these re- 
ports. 


This report begins (Section m) with a description of a hardware configura- 
tion and the software operation envisioned for an efficient implementation of an 
aided inertial navigation system. The configuration used in the flight test imple- 
mentation is also presented in this section. The runway-referenced navigation 
equations used in the experimental software and the extensions required for obtaining 
software for a global navigation system are presented in Section IV. The software 
compensation for inertial measurement unit (IMU) anomalies such as bias, mis- 
alignment and scaling of the accelerometer measurements is included in this sec- 
tion. The philosophy behind the new formulation for the error state vector for 
inertial navigation systems and its relationship to traditional procedures is pre- 
sented in Section V. The variational equations and approximations used in the 
experimental system are also summarized. The appendix gives a summary of 
vector identities and properties of orthogonal transformations which are used in 
the mathematical development. 
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H. NOTATION, DEFINITIONS AND CONSTANTS 


2. 1 Notational Conventions 

Vector and matrix notation is used extensively throughout this report. 

Scalars and vectors will be denoted by lower case symbols, and any ambiguities 
of usage between the two will be resolved in the text. Components of vectors 
will be denoted by right superscripts. Coordinate reference frames for 
vectors are indicated by left subscripts on the vectors. For example, 

12 3 

x = position vector in the "c" frame, with components x , x , x . 
c c c c 

Upper case symbols will denote matrices. In particular, the orthonormal 

3x3 matrix for coordinate transformation from the "n" frame to the "c" 

frame is denoted by T . An example of such a transformation is 
c 

n 

x = T x . 
c c n 

The notation of " • " over a symbol will have its customary meaning 
of differentiation with respect to time. The " * " (hat) mark over a symbol 
means "estimated" or "computed" value of the symbolized quantity, while 
the " ~ " (tilde) indicates that the quantity is an error or small variation. 

For instance, if x is the true value of position, it may be written as the 
sum of the estimated position and the position error. 

X = X + -x.- 
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2. 2 Reference Frames 


All reference frames used in this report are right-handed orthonormal 
frames. The following lower case letters are used for their identification, 

c computer frame, fixed to the Earth at a runway or waypoint. 

e Earth-fixed equatorial frame with "1" axis along the Earth's 

polar axis and "3" axis in the equatorial plane at the Green- 
wich meridian. 

i inertial frame, non-rotating with respect to space. 

local level reference frame at the aircraft's position. 

n navigation frame, "tied” to the platform by a coordinate 
transformation. 

p platform frame, fixed to a stable platform with "1” axis 
along the #1 accelerometer axis. 

r runway frame, fixed to the Earth with "1" axis along the 

runway and "3" axis along the local vertical at the runway. 

Figure 1 shows some geometrical relationships between the "e" frame 
and two local level reference frames. One of these " l" frames is a north- 
pointing frame and the other is a "wander azimuth” reference. The computer 
and runway frames are also level frames, but these are fixed to the Earth 
and do not change orientation with changes in the aircraft's position. 
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"1" Earth's polar axis 




Key 
1 , 2 , 3 
n, w, u 
i, j, u 
</> 

X 

a 

h 


basis vectors in Earth-fixed frame 

basis vectors in local level north -pointing frame 

basis vectors in wander azimuth frame 

geodetic latitude 

longitude 

wander azimuth 

altitude 


Figure 1. - Geometrical Relationships Between Earth-Fixed and Local Level 
Reference Frames 
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2. 3 Definition of Symbols and Constants 


Roman Symbols 


a 

b 

C 

a 

C 

g 

c.. 

i) 

d 

g 

d r* 
x 

d r* 

y 

A 

dx 

e 


f 


f 

a 

*b 

F 

n 

F 

z 

g 


h 

h 

z 

I 

M 

n 


equatorial radius of the Earth (6378. 16 km) 
polar radius of the Earth (6356. 77 km) 
accelerometer compensation matrix 
gyro compensation matrix 

elements of the transformation between local velocity and craft rate 
force-dependent gyro drift compensation 

differential operator defined by equation (5. 28) 

differential operator defined by equation (5. 29) 

incremental state estimated by the filter 

eccentricity of the oblate Earth model (. 0818349) 

T 

unit basis vector, e = ( 0 0 1 ) 

O 

(1) specific force vector (compensated delta-velocities from accelerometers) 

(2) flattening of the oblate Earth model (. 00335363) 

measured specific force (raw delta-velocity data from the accelerometers) 
accelerometer bias 

Jacobian matrix of the state rates with respect to the random forcing functions 
Jacobian matrix of the state rates with respect to the estimated state 
gravitational acceleration 

/ 2 

gravitational acceleration at the equator, zero altitude (. 00978027 km/sec ) 
altitude (usually above mean sea level) 

first-order altitude approximation from "c" frame position coordinates 
3x3 identity matrix 

Jacobian matrix of Earth central rotation with respect to position 
white noise vector (random forcing function) 
white noise on the accelerometer bias model 
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n white noise on the gyro drift model 

r^ apparent radius of curvature of the Earth 

r meridianal radius of curvature of the Earth 

m 

r normal radius of curvature of the Earth 

n 

ft rotation matrix 

t time 

T matrix of transformation between coordinate reference frames 

t. time at the beginning of a Kalman filter cycle 

t.. elements of the transformation matrix, T 

ij 

v velocity of the aircraft 

x (1) position of the aircraft 

(2) state vector of the Kalman filter 

x^ east component of aircraft position 

x level component of aircraft position 

x^ north component of aircraft position 

y measurement vector from the navigation aids 

z error state vector 

Greek Symbols 

a wander azimuth angle, platform azimuth 

a.^ azimuth of the runway or RNAV airway 

j3 rotation vector of tilts between the "c" and "n" frames 

y Earth central angle between the "c" and " V' frame origins 

r Jacobian matrix of gravity with respect to position 

0 rotation vector for Earth central angle between "c" and "l" frames 

0 unit rotation vector, 0 

u 

X longitude 

v rotation angle in the local level plane 
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p rotation vector 

Tj time constant for the correlated noise model of accelerometer bias 

T, duration of one Kalman cycle 

K 

T time constant for the correlated noise model of gyro drift 

CO 

<p rotation vector of platform tilts 

$ state transition matrix 

state sensitivity to random forcing functions 
ip geodetic latitude 

jj) geodetic latitude of the runway or waypoint reference center 

(0 rotation rate of reference frame with respect to the Earth 

craft rate, desired platform rate 

(o, drift rate between the "c” and ’’n” frames 
d 

co^ sidereal rate of rotation of the Earth (15. 041067 deg/hr) 

co platform command rate compensated for gyro scaling and misalignment 

& 

CO steady gyro drift rate 

to leveling control rate of the platform 

CO maximum leveling control rate 

max 

total rate of rotation of the platform 


\ 
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m. 


AIDED INERTIAL NAVIGATION SYSTEMS 


3. 1 Overview 

Figure 2 is a block diagram of the aided inertial navigation system 

configuration considered in this report. The main parts of the system are: 

* 

(1) the inertial measurement unit (IMU) which consists of a stable platform 
on which three integrating accelerometers sense the specific force vector 
acting on the case in the form of "delta-velocity" increments (actually, integrated 
specific force over a fixed time interval); (2) NAVAID receivers and transducers 
which provide the external "aiding" data; and (3) a digital computer, containing 
logic for IMU compensation, the navigation equations, and a modified Kalman 
filter for processing the external data. 

The raw accelerometer data is corrected for biases, alignment errors, 
and scale factor calibrations in the IMU compensation section of the computer 
logic. This compensated data is used with the computed gravity vector and 
coriolis force to obtain the aircraft's acceleration which is then integrated in 
the navigation equation logic to form current position and velocity estimates. 

The navigation equation logic also calculates the desired platform leveling 
commands, which are modified in the IMU compensation logic to account 
for gyro torquer scale factors, gyro axis misalignment, and gyro drifts. The 
resulting signal feeds the gyros which, in turn, cause the platform drive 
systems to follow. 


* Aided inertial navigation systems can incorporate strapped-down platforms too, 
but the discussion is restricted to stable platform implementations in this report. 
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INERTIAL MEASUREMENT UNIT DIGITAL COMPUTER 

(IMU) 



Key: 


f a = digitized accelerometer data 

c f = compensated accelerometer data in computer frame 
w c = desired platform rate 
^ = compensated platform command rate 
x(t) = estimated state 

dx(t) = incremental state estimated by the filter 
y(tj)= vector of measurements from the navigation aids 


Figure 2.- Block Diagram of an Aided Inertial Navigation System 








Figure 3 illustrates a sequence of computer calculations suitable for 
an aided INS such as that considered in this report. At discrete times, 
measurements from the navigation aids (e.g. , barometric altimeter, VOR/DME, 
and ILS receivers, etc. ) are keyed into the Kalman filter section of the computer. 

The current state estimate, x , is also keyed into the filter at the same time. 

The Kalman filter forms measurement residuals by computing estimates of each 
of the measured quantities based on x and subtracting these from the actual 
measurements. The filter calculates an incremental state, dx, which is an 
estimate of the error in the state estimate, x, based on the residuals, the 
measurement error model and the INS error model. The elements (components) 
of the vector, dx, typically include position, velocity, "tilt" angles, and perhaps 
an assortment of biases (e. g. , accelerometer bias, gyro drift, measurement 
biases). At the end of the Kalman filter cycle, the computed incremental state 
vector is added to the estimated state vector in the navigation equations, the 
measurement model, and the IMU compensation logic. The accelerometer data 
is processed at a high rate, so the estimated state is nearly continuous. The 
computer time required for the filter calculations depends on the complexity 
of the error model used in the filter and upon the amount of external data to be 
processed by the filter. A data rate much lower than that used for t£ie accelerometer 
data is usually more than adequate for effective aiding, and the filter need be 
cycled at only a relatively modest rate. 

3.2 The RAINPAL System 

The RAINPAL configuration is like that shown in Figure 2, except that the 
INS, a Litton LTN-51, is used strictly as a source of accelerometer data. 

A separate computer (an SDS 920) performs all the aided inertial navigation 
Computations. The SDS 920 computer estimates platform tilts and takes them 
into account in the navigation calculations, but does not torque the platform. 

This torquing is commanded by the LTN-51 computer on its own, operating as 
an independent free inertial system. The following data is transferred from 


- 11 - 



( 1 ) 


( 2 ) 


(3)(4) (5) ( 6 ) (7) ( 8 ) 


1 I 1 1 1 

L_\J 

1 i 1 'st 

rlWWWwNHIIII 

III! 

5 T k »4-. . 5 T 

1 1 1 II 1 1 1 1 1 

Niii’ 


Solve navigation equations using accelerometer 
delta velocity measurements 


Events 

( 1 ) Save state estimate and measurements 

(2) Compute incremental state and update square-root 

covariance for the measurements at time t. + r. 

i k 

(3) Save state estimate at the midpoint of the filter cycle 

(4) Compute the state transition matrix for one filter cycle 
using the state estimate at the midpoint 

(5) Propagate the incremental state and the square-root 
covariance to the end of the filter cycle time C+ 7 ^ 

( 6 ) Reduce the square -root covariance to triangular form 

(7) At time t = t^ + T^add the incremental state to the state estimate 

( 8 ) Start next filter cycle 


Figure 3.- Sequence of Computer Calculations for One Kalman Filter Cycle 
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the LTN-51 to the SDS 920 computer: 

1. digitized accelerometer data at the repetition rate of 20 Hz, 

2. aircraft attitude with respect to the platform on request, 

3. serial digital data giving latitude, longitude, true heading, and 
north-south and east-west components of velocity. 

The accelerometer data is obtained with sufficient resolution to permit RAINPAL 
to navigate "free inertial" with accuracy equal to that of the LTN-51. Attitude 
data is used in RAINPAL only to apply second-order corrections to certain of 
the external aiding data, a function which does not require high precision. Of 
the serial digital data, only true heading is used in RAINPAL (for initialization), 
but the other quantities serve as an indication of the performance of the inde- 
pendent LTN-51. 

The RAINPAL software processes accelerometer data at a 10 Hz rate. 
The error model used has 11 state variables (three positions, three velocities, 
three "tilts," plus barometric altimeter and vertical accelerometer biases). 

The resulting computer time (on the SDS 920) is such that the filter can be cycled 
at a maximum frequency of . 5 Hz. 


•r A '-*v 


V 


Js&L 


r t ■ ! 
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IV. NAVIGATION EQUATIONS AND IMU COMPENSATION 


4. 1 General Considerations 

The coordinate frame selected for the navigation equations has a 
significant effect on the calculation time and numerical precision required 
by a navigation system, both in the internal calculations and in the provision 
of appropriate external information for pilot or autopilot usage. In the 
internal calculations, we have the requirement for processing accelerometer 
and navigation aid measurements and for calculating platform commands. 

The external calculations carry a requirement for providing appropriate 
displays such as: 

(a) the position of the aircraft expressed in the global coordinates 
(geodetic latitude , longitude and altitude) and in relative coordinates 
(distance to/from en route waypoints or runway threshold, cross- 
track error and altitude) and 

(b) the velocity of the aircraft with respect to the ground, expressed 
in a local reference (north-south velocity, east-west velocity and 
altitude rate). 

7 8 

Commercially available inertial navigation systems ’ generally use a 
g 

wander azimuth coordinate reference for velocity calculations. This is 
a local level reference and north-south and east-west velocities are calculated 
from wander-azimuth components of velocity by a coordinate transformation. 
The position is calculated and retained in a transformation from an Earth-fixed 
equatorial reference frame (in which one basis vector lies along the earth's 
polar axis and another is in the equator and plane of the Greenwich meridian) 
to the local level reference at the aircraft. . The latitude and longitude are 
calculated by the use of inverse tr igonorinBi^pic tions from selected elements 
of the transformation. The altitude and afuu5ae*rate are not generally obtained 
from the INS. The commercially available inertial navigation systems have a 
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relative position output capability in the area navigation (RNAV) reference 
frame consisting of distance to/from a waypoint, cross-track error and 
altitude. Great circle paths between waypoints define the desired RNAV 
airway. 

The wander azimuth reference frame has several disadvantages when 
we consider its usage in an aided inertial navigation system for all phases 
of flight. These are: 

(1) The position is normally retained in global (rather than relative) 
coordinates. The high numerical accuracy required during the 
landing phase necessitates added calculations in fixed-point computers, 
such as double-precision arithmetic or automatic scaling. 

(2) The internally carried components of velocity (and position) require 
numerous calculations before they are appropriate for pilot display 
or autopilot usage. 

(3) Processing of navigation aid measurements is not as efficient as 
desired, since these measurements are usually relative-type 
quantities such as range and bearing to VORTAC stations. 

A convenient reference frame for navigation during landing is one which 
is aligned with the runway and centered at the threshold. In this frame, the 
internally calculated variables are distance to threshold, cross-track error, 
altitude above threshold and their corresponding velocities. The selection of this 
reference frame therefore simplifies the external calculations, minimizing the 
calculation time and numerical errors during the approach and landing phases of 
flight where rapid and accurate navigation is required. 

If we were to select an Earth-fixed reference frame which is aligned 
with the desired RNAV airway and centered at the waypoint, then distance 
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to/from the waypoint, cross-track error and a vertical distance would be 
directly available along with their corresponding velocities. The latitude 
and longitude can also be output with computation by series approximations 
which are potentially simpler than the inverse trigonometric functions used 
in the wander-azimuth mechanization. A further advantage in selecting this 
computational reference frame over the wander-azimuth mechanization is 
that processing VGR/DME data from stations in the vicinity of the waypoint 
involves simpler calculations and higher numerical accuracy with only 
single-precision arithmetic. 

The runway (or waypoint) type of coordinate frame is adopted for the 
calculations in the system described in this report. Coordinate center 
shifts which involve a translation and rotation of the coordinates as the 
aircraft proceeds from one waypoint reference to the next will be described 
in the development. 


4. 2 Navigation Equations 

The vector equation for acceleration of the aircraft with respect to the 

Earth (in a reference frame which is rotating at co + o> with respect to inertial 

, . 9,10 4 * 6 

space) is 


v = f + g - ( 2 + w) x v 


(4.1) 


where 


u> = 

u) = 
e 

f = 
g = 


rotation rate of the reference frame with respect to the Earth 
rotation rate of the Earth with respect to inertiai space 
specific force 

gravity, including centripetal acceleration. 


- 16 - 



If we select a local tangent plane frame at the take-off runway to start and 

at discrete times use a center-shift plus a rotation of coordinates, then co = 0. 

The transformation which transforms vectors from the aircraft's local level 

frame to the computer "c" frame (runway or waypoint) is T^: The gravity 

* c 

vector is represented very simply in the "t," frame. We may rewrite (4. 1) 
more specifically by including the reference frame in which each vector is 
expressed (denoted by a left subscript). 


v 

c 


T \ T P f 
c X p 


+ S) 
l 


2 oo x v 
c e c 


The position vector with respect to the coordinate center satisfies 


(4.2) 


x = v. (4. 3) 

c c 

The calculation of the quantities in (4.2) is simple except for the transformation, 

l 

T . The specific force vector, f, is obtained from the accelerometer measure- 
C P 

ments as shown on Figure 2. The platform-to-level transformation, ^T p , is part 

of the IMU compensation model covered in Section 4. 4. The gravity vector, g , 

can be approximated by a scalar function of altitude and latitude which acts along 

the vertical direction in the " V' frame. The Earth rate vector, to , has a con- 

c e 

stant magnitude of about 15 degrees per hour and lies along the Earth's polar 
axis direction. This vector is a constant in the "c" frame and changes only at co- 
ordinate center shifts. To obtain the T^ transformation, we note that what is 

c 

desired is a transformation from a level frame at the aircraft to a level frame at 

the coordinate center. Neglecting azimuth considerations for now, we can calculate 
l 

T by a single rotation in the plane formed by the local vertical vectors of the "c" 


* One of the state variables in the present RATNPAL system is the bias of the 
vertical accelerometer. The barometric altimeter measurements calibrate 
this bias continuously and, as a result, g variations have only insignificant 
effects on the accuracy of the system. 
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and " V ’ frames. Calculation of is complicated by the fact that the Earth is 

c 

not spherical. An approximate solution is to form a unit rotation vector, 0^, 
defined by 

0 = x x e / 1 x x el (4. 4) 

u c 3 lc 3 1 v ’ 

T 

where e = (0 0 1 ) and x is the position vector in the computer frame 

o C 

(i.e. , the aircraft's position with respect to the "c" frame origin). Let y be 
the rotation angle (to be defined) and then define the vector , 0 , by 0 = 0^siny . 
Then the transformation from level to computer frame is 

= I - (0x) + (0x ) (0x )/(l + cosy ) . (4.5) 

We approach the problem of defining the rotation angle, y, by allowing two 

radii of curvature for the Earth to approximate the oblate spheroid. One is for 

rotations about an axis which is parallel to the east-west direction at the center 

of the reference and the other is for rotations about an axis parallel to the north- 

south axis. The two radii are shown as functions of latitude on Figure 4. The 

formulation of the equations is given in Reference 9. The difference in radii is 

largest at zero latitude, whereas the change with latitude is largest at 45°. The 

approximate radii of curvature may be obtained by a Taylor' s series expansion 

of the exact expressions (see Figure 4), 

2 

r = r (0) + (3 e sin>J) cos ib ) (x /2) (4.6) 

m m r r r n 

and 

2 

r = r (0) + ( e sin ; b cos il) )(x 2) r (4.7) 

n n r r r r n 

where 

1 2 

x^ = north component of aircraft position ( ^x cosc^ - ^x sino! r , 

where is the azimuth of the runway or RNAV airway), 

e = eccentricity of the oblate Earth model, 

= latitude of the runway or waypoint , 
r 

r (0),r (0) = radii of curvature at the runway or waypoint, 
n m 
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radii, km 


6400 



2 9 

r = meridianal radius of curvature = a(l - e*)/ (1 - e sin^ i/> ) 


3/2 


2 2 1 

r = normal radius of curvature = a/( 1 - e sin )i> ) 
n 


Figure 4.- Radii of Curvature Versus Latitude 
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The apparent radius of curvature, r , is calculated by the (empirically formulated) 

a 

expression 


r 

a 


1 / 


2 . 2 2 
x /( x r ) 
l " 


n 


n 


2.2 
x /(x r 
e t 


m 


(4.8) 


where 

2 1 2 2 2 

x. = ( x ) + ( x ) , square of the horizontal distance, and 
-•C c c 

2 2 2 

x = x - x , square of the eastward position component. 
g n 

The altitude may then be found using 


h = h (1 - 0. 5h/r ) 
z z a 


(4.9) 


where 

3 T 

h = x+ x x / (2r ) , the first-order altitude, 
z c c c a 

The sine of the rotation angle is found using 

sin y = x / (r + h) . (4.10) 

a 

A numerical check of the approximations given in (4. 9) and (4. 10) shows errors 
smaller than one meter and 10 microradians, respectively, for distances less 
than 500 km from the coordinate center and altitudes below 20 km. The errors 
vary with direction of the aircraft from the runway or waypoint. 

4. 3 Platform Leveling Commands 

The effects of many platform error sources can be minimized by keeping 
the platform in a near-level condition. The rotation rate (with respect to inertial 
space) required for accomplishing this is 

°°l = °e + % ( 4 - 1]L) 

where 

- rotation rate of the Eartli with respect to inertial space, and 
W — rotation rate of the craft with respect tc the Earth. 
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In a local level reference frame the craft rate, , has two components in the 
level plane. The rotation rate about the vertical may be specified at any value 
for convenience. Normally it is taken as zero and the resulting mechanization is 
called a "wander-azimuth" mechanization, since the azimuth of the platform at 
any time depends on the path the aircraft took to arrive at its current location. 

For reasons to be discussed later, we will specify the component about the vertical 
such that the azimuth rate about the vertical in the runway (waypoint) frame is zero. 


In the local level frame, the craft rate is related to velocity as follows. 


1 ■ 
f°c 


C 11 

°12 


l “ 

l v 

U C 


_ C 21 

<M 

O 


2 

. l V 


(4. 12) 


The elements, c.. , are define by 
i] 


C 11 “ C 22 “ 2f ^ t ll )( f, t 21 /a 


C 12 + 2f ^W -h / a ^/ a 

2 2 

= [l-f(.t®) + 2f(,t® ) - h/a ] / a 


where f is the flattening of the oblate Earth model and a is the equatorial radius. 


The elements of the ^T 


transformation have been denoted 


,t?.. 

t lj 


T 6 = T C T 6 
l l c 


(4. 13) 


These elements can also be written as functions of latitude, 0 , and platform 
azimuth, a . 

A e i 

t = cos ! p cos a 

*0 J ..1 

^t 21 = -cos ip sin o' 

t‘31 = Sin * 
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A derivation of the coefficients, c.. , and a discussion of the validity of the ap- 
proximations are contained in Reference 10. 


The third component of ^ 
respect to the "c" reference is zero. 


is defined such that the azimuth rate with 
If we. let [see Eq. (A. 40)] 


R 2 R 3 ] 


“ - [ »l'*l + « 2 R l“2 + 9 3 R 1 R 2“3 ll! i TC < 4 ' 14) 

3 

then i> c is selected to cause 0^ (the azimuth rate) to vanish. Equation 

(4. 14) can be used to select this value of . 

I c 


3 

co 

c 


a 2 p c . c -» 2 

tan = ^‘ssVssVc 


(4. 15) 


2 

This quantity is very small compared to ,u> so long as the "c" reference 

'v C 

is within 500 km of the craft (at which range |tan 0^ | s . 0786). This particular 
platform command law will be called an "Earth-fixed azimuth" mechanization. 

It has the property that the azimuth on landing is the same as on take-off regard- 
less of route, speed, and time, if the "c” frame is held at the same reference 
point throughout the flight. 


The JT transformation (equation 4. 5) has no rotation about the vertical axis 
(azimuth). It is necessary to maintain the proper azimuth in transforming the ac- 
celeration measurements from the platform reference to the runway or waypoint 
reference. As a result of the torquing about the vertical axis, the azimuth problem 

in the transformation, T^ = T ' .T* 5 , is resolved by a fixed azimuth in .T^. This 

c cl J l 

point may be clarified by considering the sketch, 

^ ^ ' (3) 



coordinate 

center 



i 


y 


craft 


(1) 
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With either the "Earth-fixed azimuth" or the "wander azimuth" mechanization 
the platform stays fixed with respect to the Earth when the craft is stationary. 
When the craft is moving, as shown in the sketch, there is no torque about the 
vertical axis in the level frame with the wander azimuth mechanization. As a 
result, the platform's orientation with respect to an Earth-fixed frame will 
"wander" and take on a value which is dependent on the path. With the "Earth- 
fixed azimuth" mechanization, the platform is torqued about the vertical in the 
level reference to maintain a zero rate about the vertical at the center of the 
Earth-fixed reference frame. The net result is that a simple rotation transfor- 
mation transforms level to Earth-fixed coordinates. A constant rotation about 
the vertical transforms the platform reference to level axes (except for tilts). 


4. 4 IMU Compensation Model 

A simple block diagram of the IMU compensation model is presented in 

Figure 5. The accelerometer outputs are compensated for known (or estimated 

by the Kalman filter ) biases, scaling, and misalignment. The three-vector, f , 

P 

of Figure 5 is the compensated specific force* in the platform reference frame. 
The tilt transformation, T^ , transforms the specific force vector to the local 
level reference frame. The tilt transformation obeys a differential equation which 
is updated in the box labeled "Tilt Update." The leveling control rate, to^, 
drives the platform to a level condition when tilts are detected. This control law 
will be discussed in the next section. The IMU command rate is a compensated 
signal which drives the platform to a local level orientation and holds it there. 

The compensation includes scaling and misalignments (between platform and gyro 
input axes), compensation for steady drifts and drifts resulting from mass 


* Accelerometer outputs are in units of velocity since the digitizing circuitry is 
a voltage/frequency converter and a counter. By dividing the output by the time 
interval of counting (a constant), the result is an average specific force. Since 
the time intervals are short (. 1 sec or less), one may conceptually view the out- 
puts as indications of instantaneous specific force. 
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Compensation Compensation for 

for Bias Scaling and Compensation for Compensated 



Figure 5. - Block Diagram of Software Compensation for IMU Anomalies 









unbalance and anisoelastic effects. These latter drifts are functions of the specific 
force and may be compensated by an appropriate model. 

Figure 5 shows a fair degree of sophistication in the software compensation 

of IMU anomalies. The amount needed in any specific situation is obviously dependent 

on the size of the hardware anomalies and the accuracy requirements of the application. 

5 6 

In the flight test validation of the RAINPAL ’ navigator, an LTN-51 navigation system 

was used for accelerometer outputs as was mentioned earlier. The LTN-51 operated 

in its normal free-inertial mode and separate outputs of the platform accelerometers 

fed the RAINPAL software. Since the RAINPAL software could not drive the LTN-51 

platform, no calculations of command rates were necessary. The accelerometers 

t 

were compensated for biases and the tilt transformation, T , was used to com- 

^ 

pensate for the platform tilts. The transformation update was omitted since T is 

P 

zero when one cannot torque the platform. The results obtained in these flight tests 
demonstrate the excellent performance of the LTN-51 INS hardware when used as part 
of an aided inertial navigation system. 

4. 5 Platform Leveling Control Logic 

The software for compensation of IMU anomalies shown on Figure 5 includes 

l 

a compensation for platform tilt with the transformation, , such that it is not 

necessary to maintain the platform in an exact level condition. A leveling control 

l 

rate, U) , is used to drive the platform and the T transformation to level as 

4/ P 

"tilts" are estimated by the filter. The control logic is entered at the beginning of 
each Kalman filter cycle. This logic implements the following steps. 

1. Calculate the rotation vector, p, of a rotation matrix, R, (see equation A. 44) 

l 

such that R T represents a transformation involving only a rotation about 
P 

the vertical axis: 

“1 £ 

[i - (px) + (px)(px)/( 1 + cos sin | p|)]^T 
or 

R(p, sin _1 |p|) T 4 = R(e 3 ,.v) 


cos V sin v 0 

- sin v cos v 0 

0 0 1 . 

(4. 16) 
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The right-hand side of equation (4. 16) is a rotation about the local 
vertical. The platform is not commanded with leveling torques on this 
axis. 

2. ■ Calculate the control rate, u) > which will drive the platform to level 
in one filter cycle (period = r, ). 

K 

0J t = p/T k (4.17) 


Check whether the control rate is too large on any axis according to 


if to' 


to 


max 


for i = 1 , 2 


then set u? = /J * 03 / 1 to 1 for j = 1 , 2 

t. % max I j J 


(4. 18) 

(4. 19) 


Only two components of to. need be checked, since the third is zero as 

*0 

a result of the control law of step (1) as will be shown, 

4. Insert to in appropriate cells for controlling the platform and the 
» ^ 

T transformation. 

P 

In step (1) we see that a value of the rotation vector, p , is to be found to 
make (4. 16) hold. An examination of the right hand side of (4. 16) shows that 
the third row of R ( p , sin _1 |p|) must be equal to the third column of T^. The 
orthonormal properties of the matrices will then cause zeroes and unity to 
occur in the proper locations in the product. This constraint and the expansion 
of R( p , sin -1 p ) yields the following equations. 


2 

P + 

pV/( 1 

+ cos sin p |) 

= t^ 
p 13 

(4. 20) 

1 

-p + 

P 2 P 3 /(1 

+ cos sin *| p |) 

= t^ 
p 23 

(4.21) 

1 - 

[ (p 1 ) 2 + 

2 2 1 

(p ) J/ (1 + cos 

sin 'Vl> " p‘33 

(4. 22) 
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The solution of the simultaneous equations, (4.20) and (4.21) with p =0 yields 


- t^ 
p 23 


v' 

p 13 
0 


(4. 23) 


which is the solution of step (1 ) of the leveling logic for the rotation vector. 


4. 6 Coordinate Center Shifts 

For en route applications the navigation system is initialized relative to 
the take-off runway coordinate system, waypoint #1 of the sketch. 



Waypoints 2 through 6 define the en route RNAV path. The system switches to 
new waypoint-centered reference frames as the flight progresses. The spacing 
between waypoints must be less than 500 km for good free inertial performance. 
The system will insert waypoints as necessary to keep the coordinate centers 
within the limits. Waypoint #7 is the coordinate center at the landing runway. 

The system switches coordinate centers when the next waypoint in the se- 
quence is closer than the current waypoint. Prior to switching, it is assumed 
that the transformations relating waypoint "i" with "i+1" (i. e. , the new "c" 
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frame), - + ^T l , and the coordinates of waypoint #i+l in the "i"-centered reference, 
.x , are available. Then the equations for shifting position and velocity between 
centers are 


i+1 


- . ,T ( .x - .x. \ 

1+1 v 1 1 1+1' 


and 


i+1 


= . ,T C.v) 

i+1 ' l ' 


(4. 24) 
(4. 25) 


A new transformation is required for in order to transform the 

accelerometer measurements to the new center. This is found from 


tpP _ rpi ipP 
i+1 i+1 i 


(4. 26) 


The transformation, T^ , is calculated using . x in the transformation 
1+1 1+1 

described in section 4. 2. A new transformation for T 1 * is then calculated 

l 


such that 


p i+1 p 

T p = T . ,r . (4.27) 

l l 1+1 


The coordinate center shift is completed by providing Earth rate in the "i+1" 
frame. 
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V. VARIATIONAL EQUATIONS 


5. 1 General Considerations 

11 12 

The traditional procedure ’ used in the development of linear 
variational (error) equations is to define an "ideal" INS as one which 
incorporates perfect components and a perfect set of navigation equations. 
When such a system is given the true initial conditions, it will maintain the 
true state (position and velocity) for an indefinite time period. For stable 
platform systems, the ideal INS includes a perfect platform which is kept 
level at the true location by the navigation equations. The dynamic equations 
(mathematical model) for the ideal system are formulated in a fashion similar 
to those of the actual error-prone system. The error equations are formed 
by subtracting the actual equations from the ideal equations and setting the 
result equal to the mathematical model of the various hardware anomalies 
and software approximations which cause the difference. In those situations 
where the model is nonlinear, a Taylor's Series expansion is performed and 
only the first order (linear) effects are retained in the error equations. 

The primary use of the ideal INS in the above development is that it 
is an "errorless" navigation system. Another definition of an "errorless" 

INS is one where the actual hardware contains anomalies which are perfectly 
compensated by a perfect "inverse" model of the anomaly in the onboard 
computer. With reference to Figure 5, we have shown compensation for such 
anomalies as (1) bias, scaling and misalignment of accelerometers, (2) steady 
and specific-force-induced gyro drifts, and scaling and misalignment of the 
gyros, and (3) tilts of the platform. If this IMU compensation model were 
perfect and a perfect set of navigation equations were implemented with the 
true state in the computer, then we would have another "errorless" navigator. 


-2.9- 


It remains errorless even when the actual platform is tilted away from the 
vertical because these tilts (as well as the other platform anomalies) are 
exactly known and properly taken into account in the computer. With this 
philosophy of an "errorless” navigator, one tends to blame actual system 
errors on imperfect computer compensation rather than imperfect hardware. 
For example, platform tilts are not errors in themselves; rather, the error 
in the INS arises because of the difference between the true tilt and the estimate 
of tilt which is used in the computer's compensation. 

If we use the above definition of the "errorless" INS in place of the 
"ideal" INS, v/e can develop the error equations in the manner described 
previously for the "ideal" INS. If we retain the same coordinate frames for 
the error equations in the transition between philosophies, the resulting 
error equations are the same. This "errorless" concept has the advantage 
of leading to a more straightforward derivation of the error equations, as 
will be seen in the sequel. Another advantage of this approach is that it lends 
itself more naturally to aided INS applications in which tilts are estimated 
using external data. This may be seen by referring to Figures 2 and 5. If 
the filter state vector (dx of Figure 2) contains IMU errors such as vertical 
accelerometer bias and tilts, then incremental estimates from the filter can 
be used to instantaneously change the computer compensation model of these 
IMU anomalies. The effects of the error in the actual navigation equations 
are generally reduced as external data is used to estimate such anomalies 
and the behavior of the actual system tends to approach that of the perfectly 
compensated "errorless" INS. 

The error equations developed by either of these approaches contain 
forcing functions involving a product of tilt error and the measured specific 
force in the differential equation for the velocity error. These functions create 
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some problems in implementing aided inertial systems since the measured 
specific force is not well behaved and can have rapid changes in accelerometer 
outputs caused by noise and vibration. The forcing functions arise in the 
error equations because of differences between the coordinate frame of the 
measurement and the coordinate frame for the computer calculation of the 
craft's acceleration (and velocity). Traditionally the "errorless" or "ideal" 

INS takes the computer reference frame as being correct (for the development 
of error equations) and blames the platform for these differences. We introduce 
instead the concept that the craft's acceleration and velocity should be calculated 
in a reference frame which is mathematically tied to the platform. The 
result is that the specific force is in the proper frame and tilt errors between 
the platform and computer frames do not exist. Instead, errors arise in 
the calculated acceleration because the gravitational attraction and coriolis 
force are not resolved into the proper reference frame. Errors also arise 
in the rate of change of position because the velocity reference frame is 
mathematically tied to the platform and the position reference frame is defined 
with respect to the Earth. Errors in knowledge of the velocity reference frame 
therefore cause position errors since the computed velocity cannot be precisely 
resolved into the position reference. 

The development of the error equations in the sequel will use an 
"errorless" INS with the perfectly compensated IMU anomalies and perfect 
navigation equations. The navigation equations express the acceleration and 
velocity in the reference frame which is tied to the actual platform. In the 
errorless equations, the platform orientation is known and the velocity is 
transformed to the true position reference frame. Also the gravity and 
coriolis forces are transformed to the true reference frame for the acceleration. 
This new formulation has the advantage of removing the specific force from the 
error equations as was mentioned. Another advantage, as will be seen, is that 
the new coordinates for the velocity error and a new definition of "tilt" errors 
combine to render a simpler set of error equations than would otherwise result. 
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Equations for the Errorless Navigator 


The errorless navigator's differential equation for velocity (i. e. , the velocity- 

rate equation) is written in the "n" frame, which is mathematically tied to the real 

platform frame. The transformation, T , connecting the platform ("p") frame and 

n 

the "n" frame is a mathematical reality existing in the computer and is, by definition, 
errorless. The fact that we might introduce errors by trying to relate the "n" frame 
to the computer ("c") or local level ("-L") frames does not prevent us from properly 
writing errorless or actual equations in that frame for the purpose of deriving error 
equations. It is convenient to use an Earth-fixed frame for the computer frame, "c". 
If we define the q T transformation as that which the actual navigation computer 
would use to transform from platform to computer frames, we can see that in the 
errorless situation, the "n" and "c" frames would be coincident. The errorless 
navigation equations are: 


v 

n 


~p „ c -t- c 

T p f + T T ,g - (2 T 
n p n c t n 


CO 

c e 


+ 


00,) X 
n d 


v 

n 


(5.1) 


and 



(5.2) 


00, is the angular velocity with which the "n" frame drifts with respect to the "c" 
d 

frame and, therefore, ^T obeys the differential equation, 



- CO , x T 
n d n 


(5.3) 


A a n 

We will now explain the " " (hat) notation which was introduced with Tr . When 
the hat appears over a variable in the remainder of this section, that variable is 
the one used in the actual navigation system's computer to represent the corres- 
ponding errorless variable. The computer variable is generally in error because 
it is computed from "best estimate" (but imperfect) navigation information. 
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A p 

The ^T transformation is defined in the computer as the platform-to- 
computer-frame transformation, which may be formed as the product of the 
platform-to-level and level -to-computer (Earth-fixed) transformations. 


T P = T l T P 
n cl 


(5.4) 


Note that the matrix product in (5. 4) is errorless while the individual matrices 

are not necessarily so. That is, ^T^ = n T P , since the "n" frame is defined by 

means of the transformation. The individual matrices in the product may contain 

errors, because T^ cannot perfectly transform vectors from the real platform 

to the real level frame and T cannot be perfectly computed if position errors 

c 

exist. Equation (5. 4) implies 

T C = I (the identity matrix) (5.5) 


which, indeed, is the relationship used in the actual navigator. Now, to develop 

CO, , we let 
d 


T C = T P T C = T l T P T l T° 
n n p c l p i 


(5.6) 


and differentiate with respect to time to put 



in correspondence with (5. 3). 


* . 

‘c A l ~ p c * l a p c *p 'i c *p i *c 

rp = fji (jpr rp ^ fp rpr rp ipr T + p rp 

n c £ p c -Cp npi npi 


CO x T° + (0, x T° - W, x T C + CO x T C 

ncn n-on ntn nen 


(5.7) 


From (5. 3) and (5. 7), 

“d ' - ( “c 4 “t + “e - “t> (5 - S) 

where 
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CO = the computed craft rate, equations (4. 12) and (4. 14) 
c 

A 

CO = the rate used for leveling the platform (see Figure 5) 

■v 

CO = the Earth's sidereal rate of rotation 
e 

co^ = the platform's total rotation rate 


The platform's total rate can be written in terms of rates commanded by the com- 
puter and errors in the model of the gyros' torquers (see Figure 5). It is beyond 
the scope of this report to develop CO^ as a function of all the individual gyro error 
sources. With reference to Figure 5, we will assume that 

1. calibration of the gyro drifts due to mass unbalance and anisoelastic effects, 

2. alignment of gyro input axes, and 

3. scaling of the gyro torquers 


has made these gyro error sources negligible, 
will be modeled here is a steady gyro drift rate, 
axes (see Figure 5). 


The only gyro error source which 

co , referred to the gyro input 
gs 


oo 

P 


t 




CO + „C0 ) 
c e l c' 


-1 

C co + 
g gs 



Combining (5. 8) with (5. 9), we obtain another expression for CO^ . 


CO, 
n d 


-T p r w 

n L p e 



CO + C _1 
c e g 



(5.9) 


(5.10) 


C A ^ 

We may notice that if the computed .T and T transformations were errorless, 

“v P 

the only contributor to a drift between the "c" and "n" frames would be the gyro 
drifts. The magnitude of co^ will, therefore, generally be very small. If the "n" 
frame and "c" frame are closely aligned initially, they will stay that way over long 
time intervals. 


5. 3 Errors in the Actual Navigator 

The differential equations implemented by the actual navigator are of the 
same form as those for the errorless navigator. 
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(5.11) 


v 

n 


T f f 


n 


T T“ g 


• n 


x = 


n 


(2 T C to ) x v 
' n c e n 


(5. 12) 


T n = I (5-13) 

c 


As was noted before, the " A " notation symbolizes that the quantity to which it 
is applied is the actual (i. e. , estimated) value of the quantity residing in the 

* p 

navigation system's computer. The ^T is errorless by definition, while the 
other "hatted" quantities may be in error. Earth's rate is "hatless" because 
we presume to know it perfectly (or much better than the other variables) in the 
actual navigation equations. We will use the " notation to denote "error" 
and will write the following type of relationships between the errorless variables 
and their counterparts in the actual navigation equations. 


V + V 

n n 

(for vectors) 

(5. 14) 

(i - 

* 

(for rotation matrices) 

(5. 15) 


All of the errors of the actual navigator are represented in mathematical 
form as (5. 14) or (5. 15). If we assume accelerometer and gyro calibration is 
perfect except for biases and if we assume numerical approximation errors are 
negligible, the following quantities appropriately define the error state. 

vector of position errors in the "c" frame 

vector of velocity errors in the "n" frame 

vector of "tilt" errors between the "c" and "n" frames 

vector of accelerometer bias errors 

vector of steady gyro drift errors 

* (I- fix) is a first-order approximation to j^I - fix + (£x)(/?x)/( 1 + cos sin ^l/3|)l 
which is the general matrix (see appendix) for frame rotation about the rotation 
vector, 0. 


x 

c 


V = 

n 



<b = 


a) 

gs 


- 35 - 



The "tilt" error vector, 0 , as it is defined here, is the only error state which is 
not generally familiar or readily recognizable. Tilts are usually defined between 
the level and platform frames. If we denote such tilts by the vector, <p, where 

p T* = (I-Jpx)T 1 (5.16) 

we can relate j3 to <p as follows. Errors in the computer-to-level transfor- 
mation are described by a rotation error vector, 6 . 

t T° = (I- t 8x)^T C (5.17) 

By noting that 0 is defined such that 

f = (I-/x)T n (5.18) 

and using the identity 

T n T P = T l T P (5. 19) 

c n c l v ' 

it can be seen that 

"i 

J = -( c 0+c^ < 5 * 20 > 

A c 

The transformation T is calculated from the position vector in the Earth-fixed 

*v ^ 

"c" frame. Hence, ^ can be approximated as a linear combination of the position 

errors, x. 
c 

0 = M x (5.21) 

c c 

The matrix, M, which expresses the linear relationship will be defined later. 
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5.4 Linearization 


The variational equations for the navigation equations are developed by 
differencing the errorless and actual navigation equations, then linearizing the 
result by omitting second and higher order terms. Considering first the 
position equations, (5. 2) and (5.12), we develop a first-order linear differential 
equation for position error, using (5. 14) and (5. 15) to represent the error. 


x - x 
c c 


rr . n 
T v 
c n 


-\n - 
T v 
c n 


= ( I - £ x ) T n (v + v ) - T“ ; 

c c m n c n 


* n a 
v 


3x T n v 
c c n 


c* 


A n ~ 
x T v 
c n 


(5. 22) 


Bv observing that B and v are small, we drop their product. Then, using 
' c n 

(5. 13), we may write the equation in its final desired form. 


x 

c 


(I) V 



(5. 23) 


Equation (5. 23) is the vector differential equation for the position errors. It 
shows that the rate of change of position errors in the "c" frame is the velocity 
error in the "n" frame plus a term dependent on the magnitude of the estimated 
velocity and the rotational error vector, j3 , between the two frames. 


Development of the rate of change of velocity error is formally identical to 
that of the position error, so the intermediate algebraic steps will be omitted. The 
result is 


a 5 = <n T Vb - <c Sx) </♦.*> + <„*\* - < 2 c“e x V 

- [< n v*M2 c ‘V>J/ + <n }x >n a d ' 


(5. 24) 
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Equation (5. 24) is not yet in its final explicit form in terms of the error state, 

c 

since it includes the T error, 0 , the gravity error, g , and the drift rate, 

4, C *0 

u> . But 0 and „g can each be written as linear combinations of the position 
n d c 4 

error, and Wj can be expressed as a function of the tilt and gyro drift errors. 
Let us consider each of these terms. 


For consideration of errors, the transformation from computer to level 
frames can be quite accurately computed in the vicinity of the Earth-fixed com- 
puter reference center by 


,T a (I - 0x) 


(5. 25) 


where 


0 = 


r +h 
a 


- x 
c 

-1 

X 


(5. 26) 


and where r is the Earth's apparent radius of curvature and h is altitude. We 
a 

can differentiate this representation with respect to position to obtain the matrix, 
M , of equation (5. 21). 
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r +h 
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d r 

. y 

(l+d r 1 
v x 

0 


- ( l- d y r > 

d r 2 
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d r 
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, 3 
d r 
x 


(5.27) 


In (5.27), we have made the substitutions 

-1 


d r = — — r — r (r +h) 
x r + h ~ i ' a ' 

a Ox 


and 


(5. 28) 


d r 
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-2 * 

x 0 


~7T ( r „ + h ) 


r + h .•'i a 
a Ox 


(5.29) 
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12 1 2 

For practical purposes, d r , d r , d r and d r are zero and the partial 

x x 3 y y 

derivative of ( r + h) with respect to x is unity. These simplifications were 

3 

used in the RAINPAL onboard implementation. 


Reference 12 gives an evaluation of empirical, second order, and atandard 
gravity relationships. In the RAINPAL mechanization, the gravity vector is 
approximated by 


0 

0 

2 

g^( 1 - 2 h/ a + . 00526 sin i/)) 


(5. 30) 


Whatever gravity model is used may be differentiated with respect to position 
to obtain the matrix, p, of 






(5. 31) 


RAINPAL' s error model considered variations with altitude only and for this 
case, equation (5.32) holds. 


r 


0 

0 

0 


0 0 

0 0 

0 - 2g Q / a 



(5. 32) 


The coupling of vertical errors into velocity as indicated by equations (5. 24) 
and (5. 31) is small and probably could be completely omitted from any mechani- 
zation using barometric altimeter data to stabilize the vertical channel. 


The drift rate, to^ , between the "e" and "n" frames was stated In equation 
(5.10). Using (5. 18) in (5. 10) with T n =I gives 


W, = + ( to x) /3 - ( T P C~ 1 ) W 
d c e c n g gs 


(5.33) 
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In the RAINPAL equations, oo was assumed to be described by 

gs 




& 

gs 



(5. 34) 


The time constant, t , was set at 10 minutes and the magnitude of the 
white noise, n , was set to give an "rms" drift rate of about 2 deg/hr. 


Using equations (5. 21), (5. 31) and (5. 33), we can re-write the differential 
equation for velocity error, (5.24), in its desired form. 


v 

n 


(-gx)M + T*T X + [- 2 0 0 x J v 

L c c c L ce n 

+ L(-gx)- (vx)(wx)J ^ 

L c n c e ' J c 



+ 






(5. 35) 


The variational equation for 8 can be derived by combining equations 
(5. 3), (5. 13) and (A. 9) as follows. 


- ( T n - T C ) = - 2 fix = - OJ x (I - fi x ) - (I +fi x ) u>x (5.36) 

dt c n ' r d d 

0 

fi = OJ - (fix CO ,)/2 
d a 



- ( 


00 X ) fi 
c e 


A D -1 ~ 

+ ( T^C ) 00 
n g gs 


(5.37) 


The accelerometer bias was taken (in RAINPAL) as correlated noise. 


*b 




(5.38) 


WTien the vertical accelerometer bias is included as a state variable in the 
RAINPAL mechanization, it is continuously compensated by the altimeter 
measurement. 
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5. 5 Summary and Mechanization 

The variational equations may be summarized in matrix form if we define 
a vector of errors, z. 


z 


x 

V 

n„ 


eg 



(5.39) 


We may then write the variational equations of the preceding section as a single 
equation. 


F z + F n 
z n 


(5.40) 


We utilize the equations of section 5. 4 to identify the elements of F and F . 


n 


[(- gx)M+ „T*r] 


0 

0 

0 

0 

0 

0 

I 

0 


-2 w x 
■ c e 


0 

0 

0 

0 

I 


0 [ 

0 

0 


[ n ^x ] o 0 

l~ vx 

\ l Ux 

TfI r / 
o L - — i 


- to X 

c e 
0 

0 


co 


n 


co 


(5.42) 


The individual elements of F and F shown here are each 3x3 matrices. 

z n 
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The general solution of equation (5.40) can be expressed in a simple form 
if n is approximated as a constant over a Kalman filter cycle of period . 

This solution is 

z(t. + T ) = $ (t. + r. ; t. ) z (t.) + $ (t. + t. ; t. ) n (t.) (5. 43 ) 

l k l k l l n l k x ' l ' ’ 

In (5.44) 3> is the transition matrix which obeys 

* = *yt)$ with $ (t. ; t^ ) = I (5.44) 

and <i> n is the sensitivity to forcing functions which obeys 

$ = F (t) $ + F (t) with $ (t.:t.) = 0 . (5.45) 

n z n n n i i 

We assume the state-dependent elements of F (t) may be approximated over 

z 

the Kalman cycle as constants, with their value calculated from the state at 
the midpoint. Then by use of rectangular integration of (5. 44) and (5. 45) 

* (t i +T k ; V “ 1 + F z (t i +T k /2) * T k (5,46) 

and 

$ (t. + T. ; t.) a F (t. + T, /2) * T. (5. 47) 

n i k i n l k ' k 

As may be noted from (5. 41) and (5. 42) a large number of zeroes occur in the 
two matrices given by (5. 46) and (5. 47). To avoid -wasting memory and 
calculation time, only the non-zero elements and non-unity elements are stored 
and calculated. Constant arrays containing appropriate indices allow calculations 
such as a prediction of the incremental state 

Z{t i +T k ) = Z( V + F z (t i + V 2) * V Z(r i ) (5,48) 
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without requiring multiplications by zeroes or storage of the complete F 
matrix. Equation (5.43) may also be obtained by use of rectangular 
integration of (5. 40) with n = 0. In the RAINPAL implementation the 
transition matrix as such is neither calculated nor used. 



VI. 


CONCLUDING REMARKS 


The material presented in this report extends the navigation equations and 
variational equations of the RAINPAL formulation to an accurate all-weather aided 
inertial navigation system for all phases of aircraft operations. The development 
has been for stable platform -type inertial measurement units. A great deal of the 
work, however, applies directly to strapped-down IMU's as well. Some of the 
features incorporated in this preliminary design for new INS software include the 
following. 

1. Relative position-type navigation is used for all phases of 
the mission so that high scaling gives good numerical ac- 
curacy with single-precision fixed-point arithmetic. In 
contrast, existing systems use an Earth-centered reference 
with which single precision results in a large truncation 
error. 

2. Relative navigation is used with the coordinate frame and 
coordinate center selected to minimize output calculation. 

This feature can markedly reduce the real-time computer 
requirements of the navigator and therefore release com- 
puter time for other aircraft subsystems and displays. 

3. The Earth-fixed azimuth platform drive holds the platform - 
azimuth fixed with respect to the azimuth at the coordinate 
center. This feature reduces the complexity of the calcu- 
lations and removes the drifting character of the "wander 
azimuth" mechanization. 

4. The platform tilt control law and tilt compensation trans- 
formation eliminates the coupling between filtering and 
platform control. This feature should eliminate some dif- 
ficulties in those areas where nonlinear effects have caused 
convergence problems in the past. 

5. The new approach for development of the variational equations 
leads to a better separation and clearer understanding of the 
effects of INS errors. The resulting simpler set of error 
equations which do not contain the accelerometer measure- 
ment in the Jacobian matrix of acceleration with respect to 
the state is felt to be a significant advance in INS error 
modeling. 
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Due to the significance of these results, the incorporation of this improved 
software in a flight test prototype should be rapidly pursued. The potential benefits 
of a single system for the complete flight envelope is a very attractive cost-effective 
concept. 

This report has focused on the navigation equation and variational equation 
software for improvements in aided inertial navigation systems. In the RAINPAL 
study effort, about 4K of memory was sufficient for aided inertial navigation, and 
time was available for other calculations. This system used a 1963 vintage computer 
(the SDS 920) which is slow by comparison to currently available computers. With 
typical "state of the art" computers and improved RAINPAL -type software, the total 
navigation function for all phases of flight will require about 25% real-time utiliza- 
tion. The percentage real-time utilization will reduce by another order of magnitude 
with the next generation of computers. Although there is certain justification for 
further software improvements, such improvements are not the road blocks to 
achieving the fast and accurate navigation systems of the future. 
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APPENDIX 


Vector and Matrix Relationships 


A. 1 Mathematical Operations for 3-Vectors 

This appendix presents definitions of some frequently -used relationships 
between vectors and orthogonal (rotation) matrices. 

The general 3-dimensional vector, a, is considered to be a column matrix 


and its transpose a row matrix. 
r a 1 1 

I 

a = 


r 1 2 3 ~i 

° L a a a J 


(A.l) 


T T 

Thus, the dot or inner product, a-b, is written a b or b a. The outer product, 
T 

ab , is defined by 
T 


ab 



■ 11 


LI 

1, 2 

1, 3~ 


a 

j— - — * 

a b 

a b 

a b 


2 

,1 2 .3 

2,1 

2, 2 

2,3 


a 

b b b 

a b 

a b 

a b 


3 

L J 

3 1 

3,2 

3, 3 


a 


a b 

a b 

ab J 


(A. 2) 


The vector cross-product, axb, is often written as a skew-symmetric matrix, 
(a x), operating on the vector, b. 


axb = 


r 11 


~ ll 



3 

2 

a 


b 


0 

-a 

a 

2 


,2 


3 


1 

a 

X 

b 

— 

a 

0 

-a 

3 


l 3 


2 

1 


a 


b 


-a 

a 

0 


= ( a x) b (A. 3) 


The matrix, (ax), pronounced "a-cross", has the property 

(a x) T = - (ax) (A. 4) 

so that 

(a x b) T = - b T (a x) . (A. 5) 

The vector triple cross-product is conveniently expressed as a difference 
between two vectors, 

a x (b x c) = b (a- c) - c (a- b) , (A. 6) 


which leads directly to a decomposition of the matrix product, (a x) (b x), 
(a x) (b x) = ba T - (a^b) I 

and an alternative representation of the operation (a x b) x. 


T T 

(a x b) x = ba - ab 


(A. 7) 
(A. 8) 
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Notice also that 


(a x) (b x) - (bx)(ax) = (axb)x 


(A. 9) 


A unit vector, u, may be defined from a by 



u = a/y/ a a 


from which we observe that u u = 1. As a special case of the triple cross- 
product, we have 

T 

(u x) (u x) = uu - I 
and, since (u x) u = 0, 


(-l)^ n ^(u x) for n odd 


(ux) 


n 


(_l)( n /2)( j _ UU T ) for n 


(A. 10) 


(A. 11) 


(A. 12) 


even 


The matrix, (I - uu ), is of special interest since 

(I - uuV = (I - uu T ), (A. 13) 

which is a property of idempotent matrices. This matrix projects any vector, b, 
into the plane normal to u by subtracting off that component of b which lies along u. 


( I - uu ) b = b - u (u- b) 


(A. 14) 


Conversely, the matrix (uu ) operating on b removes those components of b 
which do not lie along u. 


A. 2 Rotation "Coordinate Transformation" Matrices 

Let us now consider the 3x3 rotation matrix, R. The orthonormal 
property of R is 

R T R = RR T = I (A. 15) 

-1 T 

which is equivalent to R = R . The effect of multiplying an arbitrary vector, a, 

by R is to rotate a, since the magnitude of a is unchanged by the operation. 

b = Ra (A. 16) 

T T T T 

b b = a R R a = a a (A. 17) 

Let us now consider the effect of rotation on the vector-matrix operations which 

have been described. If we use the notation, a = Ra, b = Rb, etc., it is 

r • r 
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apparent that 


T , T. 

a b = a b 

r r T T T 

a b = R (ab ) R . 
r r 


(A. 18) 


and a b = R (ab ) R . (A. 19) 

r r 

In order that the vector cross-product be preserved under orthogonal transfor- 


mation, 


we must have 


ax b = R(axb), 
r r 


(^a x) = R (a x) R 
(an important result) for then 


a x b = R (a x) RR b = R (a x) b = (axb) 
r r r 


(A. 20) 


(A. 21) 


(A. 22) 


We will now state without derivation that the general rotation matrix can 
be written in the form 

T 

R(u,0) = cos 0 1 + (1 - cos 0 ) uu - sin 0 (u x) (A. 23) 
or R(u,0) = eos0(l-uu T ) + (uu T ) - sin9(ux) . (A. 24) 

It is readily observed that this representation has the orthonormal property. 

R(u,0) R T (u,0) = I (A. 25) 

To show that any general rotation matrix may be put into the form, R(u,0), 
we will present a means of solving for u and ©from a general R. Tfie matrix. 


C = R - R , 

is always skew-symmetric, even if R is not a rotation^ since 
T 

C = - C . 

If the elements of C are c.. , we can solve for u as the vector 

ij 


1 2 2 2 ’ 

c + c + c 
32 13 21 


(A. 26) 


(A. 27) 


(A. 28) 


which follows from the definition of the skew-symmetric matrix, u x , in (A. 3). 


The rotation angle, 0 , is found from 


2 sin 0 


2 , 2 
+ c . + c, 

13 21 


(A. 29) 
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trace R = 3 cos 0 + (1 - cos 9 ) =1+2 cos 0 


so that tan 9 = >Jc + c + c / trace R . (A, 

I 04 13 4± 

The vector, u, is the eigenvector of R which corresponds to the eigenvalue of 
unity since 


R u = u. 


(A. 32) 


It may be of interest to note that if u is rotated by the orthogonal matrix, M, 


to render 


u = M u 


(A. 33) 


then R^u.0) = M R(u,0) M , (A. 34) 

this result following directly from the form of R(u,0) and the earlier result 
that (^u X ) = M (u x) M T 


The "canonical" rotations about the references bases are special cases 
of R(u,9). For example, if 

"l 

u = 0 ( 


(A. 35) 


10 0 


0 0 0 


10 0 


R(u, 0) = cos 0 0 1 0 + (1 - cos @ 0 0 0 - sin 0 0 0 -1 


0 0 1 


0 cos 9 sin 0 
0 -sin 0 cos0 


0 0 0 


0 10 


(A. 36) 


A transformation of coordinate frames is frequently made up of three sequential 
canonical rotations in a given order. It is possible to duplicate the result with 
a single rotation through 0 about u, where 0 and u are obtained from the product 
matrix by the process described earlier. In this way 0 and u may be expressed 
as functions of the canonical rotation angles. Another item of frequent interest 
is an expression for the angular velocity as a function of the Euler angles 
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• (canonical rotation angles) and their rates. Let us assume that 
• ' ; R = R ( u 3 .e 3 ) R(u 2 . 0 2 ) R^i^i) - R 3 R 2 Ri - - (A. 37) 

where u. and 0. are the canonical axis and angle of the i-th rotation. It can . 

11 

be shown that there is always a vector, co, such that the rate of change of any 
orthogonal matrix, R, may be written 

R ■= . - to x R. • • (A. 38) 

If u = 0, to = 0u and we have the special-case result, 1 

R(u,0 ) = - 0 u x R(u,0) . ‘ (A. 39) 

Applying this special result repeatedly to the. derivative of R in (A. 37) we have 


R - Wi*Wr*'Wi ■ - ■ : 
55 Av R -Wi 18 * - S MV\ , 

- AV R - R 3®2 U 2 X)R 3 R - R 3 B 2< 5 iV> S 2 B 3 R 

" - <V 3 + A R 3 R 2V XR 

The vector sum in (A. 40) is equated to the angular velocity, to , in rotated 
coordinates. 


(A. 40) 


to 


= 0 u 
3 3 


0 R u 
2 3 2 


KWl 


Its expression in non-rotated coordinates is. found from to = 


to = 


• t T 

e 3 R lV 3 


• T 

e 2 R l u 2 


Vi 


T 

R to. 
r 


(A. 41) 
(A. 42) 


Another form of the general rotation matrix, (A. 24), can be quite useful. 

If we define a vector, b, as 

b = u sin 0 . < (A. 43) 

the general rotation matrix, R(b,0 ),. can be written 

\ • R(b,0) = I - (b x) .+ (b,x) (b x) r/(l + cos 0 ) . (A. 44) 

This form explicitly shows the error in approximating R(b,0 ) by ( I - b x ) 
for small angles, when (b x) (b x) is negligible to first order. 


Another substitution leads to the "Euler Parameters. " Define a vector, q, as 
q = u sin (0/2) • ... * (A. 45) 

Then, using double-angle identities for sin 0 and cos 0 in (A. 24), we obtain 
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R(q,0) = I + 2 (q x) (q x) - 2 (q x) cos (8/2) 

= I + 2 (q x) [ (q xj - (I) cos (8/2) 3 ' (A. 46) 

By defining a fourth element of q to be cos (8/2) , the vector of Euler's 
Parameters is completed. - • «. 




q 

cos (8/2) 


(A. 47) 


Another form, then, for the rotation matrix can be written as a function of q. 

R(q) = 2 [ ( q* - 1/2) (I) + ( qq T ) - q 4 (q x) ] (A. 48) 


A. 3 Derivative of the General Rotation Matrix 

We again consider the general rotation matrix in the form 

. rii . s . 

R(u,8) = cos 8 (I) + (l-cos0)(uu ) - sin 0 (u x) (A. 49) 

where u is a unit column vector. Let the derivatives of u and 8 with respect 

to any scalar (such as time) be denoted by u and 0 , respectively. Note that 

u U = u u = 0 ; " (A. 50) 

since u is constrained to be a unit vector. 

R(u,0 ) = [ - sin 8(1- uu T ) - cos 0 (u x) J 0 + ( 1 - cos 0)(uu T + uu T ) 

• • » 

- sin Q (u x) (A. 51) 

Equation (A. 51) may be put into more useful form by post-multiplying both 

T « t 

sides by R (u,0)R(u,0) = I and re-forming the matrix product, R(u,0) R (u,0), 

which will be shown to be skew-symmetric. Consider first only the coefficient 

of 0 in this product and leave the u-terms until later. 

[ - sin 0 (I - uu r ) - cos 0 (u x) ] R T (u,0 ) = - (u x) (A. 52) 

This preliminary resfllt is the same as (A. 39). It was obtained algebraically 

using the following identities 


T T 

( I - uu ) (uu ) = 0 

(A. 53) 

T 

( I - uu ) (u x) = (u x) 

(A. 54) 

(u x) (u x) = <uu T - I ) 

(A. 55) 

T 

(u x) (uu ) = I 

(A. 56) 
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Equation (A. 55) is a special case of the vector triple cross-product, (A. 6) and (A. 7). 
Using the vector identity of (A. 8), we now consider terms in u to complete the 
derivation of R(u,0 ). Omitting a great deal of intermediate algebra, 

[ (1 - cos 0 ) (uu + uu ) - sin 0 (u x) JR (u,0 ) = (cos 0 - 1) (u x u) x - sin 0 

(A. 57) 

• r r 

Combining equations (A. 52) and (A. 57) with (A. 51) to form R(u,9 ) R (u,0) R(u,0) 
we find 

R(u,0 ) = - Wx R(u,0 ) (A. 58) 

where 

• • • 

« =0u + (1- cos 0 ) (u x u) + sin 0 u . (A. 59) 

It may be observed that u, u x u, and u are mutually orthogonal. It may also be 
of interest to consider the fact that 

u T y = 0 . (A. 60) 

Although the "dot" notation usually denotes differentiation with respect to time, the 
result given here holds for differentiation with respect to any other quantity. 
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